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ABSTRACT 

This paper describes a model which can explain the observed clumpy structures of debris 
disks. Clumps arise because after a planetary system forms its planets migrate due to angular 
momentum exchange with the remaining planetesimals. Outward migration of the outermost 
planet traps planetesimals outside its orbit into its resonances and resonant forces cause azimuthal 
structure in their distribution. The model is based on numerical simulations of planets of different 
masses, Mpi, migrating at different rates, hpi, through a dynamically cold (e < 0.01) planetesimal 
disk initially at a semimajor axis a. Trapping probabilities and the resulting azimuthal structures 
are presented for a planet's 2:1, 5:3, 3:2, and 4:3 resonances. Seven possible dynamical structures 
are identified from migrations defined by = Mpi/Mi, and 9 = cipi^Ja/Mi,. Application of 
this model to the 850/xm image of Vega's disk shows its two clumps of unequal brightness can 
be explained by the migration of a Neptune-mass planet from 40 to 65AU over 56Myr; tight 
constraints are set on possible ranges of these parameters. The clumps are caused by planetesimals 
in the 3:2 and 2:1 resonances; the asymmetry arises because of the overabundance of planetesimals 
in the 2:l(u) over the 2:1(1) resonance. The similarity of this migration to that proposed for our 
own Neptune hints that Vega's planetary system may be much more akin to the solar system 
than previously thought. Predictions are made which would substantiate this model, such as the 
orbital motion of the clumpy pattern, the location of the planet, and the presence of lower level 
clumps. 

Subject headings: celestial mechanics — circumstellar matter — planetary systems: formation — plane- 
tary systems: protoplanetary disks — stars: individual (Vega) 



1. Introduction 

The discovery of disks of dust around main se- 
quence stars showed that some grain growth must 
have occurred in these systems, since this dust 
is short-lived and so must be continually replen- 
ished (e.g., Backman & Paresce 1993); the dust is 
thought to originate in a coUisional cascade that 
starts with planetesimals a few km in size (Wy- 
att & Dent 2002). Thus the lack of warm dust 
in the inner ~ 30 AU of these systems implies a 
relative paucity of planetesimals in these regions 
(Laureijs et al. 2002). This would naturally be ex- 
plained if the planetesimals here grew large enough 
for their gravitational perturbations to clear these 



regions of any remaining debris, i.e., if they grew 
to ^ 1000 km sized planets (Kenyon & Bromley 

2002) . While the presence of such planets is still 
questionable, it is widely believed that these de- 
bris disk systems are solar system analogs, and 
that the dust that is observed is produced by the 
destruction of analog Kuiper belts (Wyatt et al. 

2003) . 

The most positive indication that debris disk 
systems harbor planets comes from the morphol- 
ogy of the dust disks. All of the disks which 
have been imaged exhibit clumps and asymmetries 
which have been interpreted as perturbations from 
a planetary system: e.g., a brightness asymmetry 
in the HR4796 disk (Telesco et al. 2000) can be 
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explained by the secular gravitational perturba- 
tions of a planet on an eccentric orbit (Wyatt et al. 
1999); clumps observed in the e Eridani, Vega and 
Fomalhaut disks (Holland ct al. 1998; Greaves et 
al. 1998; Wilner et al. 2002; Holland et al. 2003) 
may be indicative of dust trapped in mean motion 
resonance with a planet in these systems (Ozcr- 
noy et al. 2000; Wyatt & Dent 2002; Quillen & 
Thorndike 2002; Kuchner & Holman 2003); and a 
warp in the f3 Pictoris disk (Heap et al. 2000) may 
be caused by the perturbations of a planet on an 
inclined orbit (Augereau et al. 2001). 

This paper focuses on the possibility that 
clumps in debris disks are caused by dust trapped 
in planetary resonances. Simple geometrical ar- 
guments show how material that is in such reso- 
nances is not evenly distributed in azimuth about 
the star it is orbiting (e.g., Murray & Dermott 
1999; see §4), and this observation has been 
exploited by several authors to use debris disk 
clumps to infer the presence of planets. In the 
studies currently in the literature, dust is trapped 
into planetary resonances when it encounters these 
resonances while spiraling in towards the star due 
to P-R drag, much in the same way that dust 
originating in the asteroid belt becomes trapped 
in resonance with the Earth (Dermott et al. 1994), 
and dust originating in the Kuiper belt might be- 
come trapped in resonance with Neptune (Moro- 
Martin & Malhotra 2002). However, Wyatt et al. 
(1999) showed that P-R drag is not important for 
dust in the bright debris disks that are currently 
known, since these are much denser than the zodi- 
acal cloud or the putative Kuiper belt dust cloud. 
This means that debris disk dust is destroyed in 
mutual collisions much faster than the P-R drag 
timescale until it is small enough to be blown out 
of the system by radiation pressure. 

Another reason why dust may be trapped in 
planetary resonances is evident by comparison 
with the solar system. About one third of known 
Kuiper belt objects (hereafter KBOs), including 
Pluto, are trapped in 3:2 resonance with Neptune 
(Jewitt 1999). More recently KBOs have also been 
found in Neptune's other resonances, including the 
1:1, 4:3, 7:4, 2:1, and 5:2 resonances Chiang et al. 
(2003). These planetesimals are thought to have 
become trapped in these resonances when Nep- 
tune's orbit expanded in the early history of the 
solar system (Malhotra 1995). In this scenario, the 



orbital expansion arises from angular momentum 
exchange resulting from the scattering by Neptune 
of the residual planetesimal disk left over at the 
end of planet formation; the current orbital dis- 
tribution seen in the Kuiper belt can be explained 
by the migration of Neptune from 23-30 AU over a 
period of 50 Myr (Hahn & Malhotra 1999). While 
this model does not completely describe the KBO 
distribution, it turns out that some of its puzzling 
features, such as the high inclination population 
of the classical KBOs (i.e., those not trapped in 
resonance with Neptune), can now be accounted 
for with more detailed models of the migration 
(Gomes 2003). Besides, the success with which 
it explains the orbital distribution of the resonant 
KBOs means that the migration of Ncptime is al- 
most certain to have played a role in the formation 
of the Kuiper belt. Furthermore, models of the 
formation of such massive planets in the absence of 
the substantial accretion of gaseous material pre- 
dicts such migration (Fernandez & Ip 1984). Thus 
if planets did form in the inner regions of debris 
disk systems, some fraction of the parent planetes- 
imals from which the dust originates is expected to 
be trapped in resonance with the outermost planet 
of that system. 

Current studies of planet migration and the 
consequent resonant trapping deal with the migra- 
tion of Neptune and its effect on the Kuiper belt 
(Malhotra 1995; Hahn k Malhotra 1999; Zhou et 
al. 2002; Chiang & Jordan 2002). This study 
is the first in a series which aims to determine 
the effect of planet migration in extrasolar sys- 
tems. In such systems, planets of different mass 
may have formed at different distances from dif- 
ferent mass stars; they also may have migrated at 
different speeds. This paper explores the effect of 
planet migration on a planetesimal disk, while the 
details of how the planet achieves that migration 
and of the dust distribution resulting from the de- 
struction of these planetesimals are left for future 
studies in the series. 

The layout of the paper is as follows. The nu- 
merical model of planet migration is described in 
§2 and in §3 this is used to calculate the probabil- 
ity of capture of planetesimals into different res- 
onances for given migration scenarios. §4 shows 
how the geometry of resonant planetesimal orbits 
causes their distribution to be clumpy and param- 
eters defining this dumpiness are derived from the 
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numerical simulations. In §5 I show how these 
results can be applied to the structure of debris 
disks using the specific example of the structure 
observed around Vega to set constraints on the 
mass and migration history of planets in this sys- 
tem. The conclusions are given in §6. 

2. Numerical Model 

The dynamical evolution of a system comprised 
of 200 planetesimals and 1 planet was followed us- 
ing the RADAU fifteenth order integrator program 
(Everhart 1985). This integration scheme is self- 
starting in that the time steps of each sequence are 
variable and arc chosen by the integrator based 
on the results of the previous sequence; substeps 
within each sequence are taken at Gauss-Radau 
spacings. 

In these integrations all bodies arc assumed to 
orbit a star of mass M*. The planet has a mass 
of Mpi while the planetesimals are assumed to be 
massless, i.e., they do not affect the evolution of 
the planet's orbit though they are affected by it. 
At the start of the integration, the planetesimals 
are assumed to have eccentricities, e, chosen ran- 
domly from the range to 0.01. Their inclina- 
tions, /, were chosen randomly from the range 
to 0.57° = 0.01 rad, and their arguments of pe- 
riastron, Co, longitudes of ascending node, fl, and 
longitudes. A, were each chosen randomly from the 
range to 360°. The planet was assumed to have 
all of these orbital parameters set to zero at the 
start of the integration; i.e., the planet starts on 
a circular orbit in the mid-plane of the planetesi- 
mal disk. The initial semimajor axes of the plan- 
etesimals, a, and that of the planet, api, were set 
according to the integration being performed. 

Planet migration was simulated by the addition 
of a force acting in the direction of orbital motion 
of the planet. The prescription for the force used 
in this paper has as its input the variable d^ar and 
results in an acceleration of: 

V = 0.5d„arv/GAf,/a3. (1) 

Since the work done by this force results in a 
change in orbital energy (sec e.g., section 2.9 of 
Murray & Dermott 1999), this causes a change in 
the planet's semimajor axis of d = o? jGM^ = 
dyar', i-e., the force defined by i) results in a con- 
stant rate of change in the planet's semimajor axis. 



The planet maintains its circular orbit {epi = 0) 
in the midplane of the disk {Ipi = 0). 

The integrator itself has been extensively 
tested, and the addition of planet migration was 
tested by repeating simulation la of Chiang & 
Jordan (2002). The resulting final distribution of 
eccentricities, inclinations and semimajor axes of 
planetesimals was qualitatively the same as that 
depicted in their figure 3. 

3. Resonant Capture Probabilities 

Planetary resonances are locations at which 
planetesimals orbit an integer p number of times 
for every integer p + q times that the planet orbits 
the star. Kepler's third law shows that resonances 
occur at semimajor axes of: 

/ , \ 2/3 

dp+q-.p = dpi ['~^ J ■ y'^) 

If the planet is migrating at a rate dpi, these res- 
onances only remain in the vicinity of planetesi- 
mals at a given semimajor axis for a short time. 
However, while close to a planetary resonance, a 
planetesimal receives periodic kicks to its orbit 
from the gravitational perturbations of the planet. 
These can impart angular momentum to the plan- 
etesimal so that its semimajor axis increases in 
such a way that the planetesimal is always orbit- 
ing at the location of the resonance; such a plan- 
etesimal is said to be trapped in the planet's res- 
onance. The resonant forces causing this trapping 
are discussed in greater detail in §4. 

The first set of integrations was performed with 
the aim of determining how many planetesimals 
initially at a semimajor axis a (in the range 30-150 
AU) get trapped in a given resonance of a planet 
of mass Mpi (in the range 1-300 M^) that is mi- 
grating at a constant rate of dpi (in the range 0.01- 
1000 AU Myr~^), where all bodies are orbiting a 
star of mass (in the range 0.5-5 A'/©). Stud- 
ies of the migration of the solar system's planets 
showed that the most important resonances are 
the 3:2 and 2:1 resonances, and to a lesser extent 
the 4:3 and 5:3 resonances. These are studied in 
individual sections below, where a model is derived 
to estimate capture probabilities using equations 
(6)- (8) and Table 1. Readers not interested in how 
this model is derived could skip to §3.4 without 
much loss in continuity. 
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The way the integrations were performed was 
to study specific values of a, Mpi and M*, then 
determine the trapping probabihty, P, for differ- 
ent migration rates chosen with the aim of obtain- 
ing P in the range 0.1-0.9. Trapping probabihties 
were determined by starting the planet at a semi- 
major axis doi below the location of the resonance 
and allowing the integration to continue until the 
planet had migrated da^ beyond the resonance. 
The parameters doi and da2 were chosen by mon- 
itoring the orbital elements of the planetesimals 
at 10 intervals throughout the integrations. Since 
resonances have a finite libration width, doi was 
chosen so that the planetesimals were started well 
outside the resonance region. At the end of all 
runs, da2 was chosen so that two distinct pop- 
ulations of trapped and non trapped planetesi- 
mals were easy to distinguish by the semimajor 
axes and eccentricities of their members. Typ- 
ically both doi and da2 were set between 0.5-4 
AU. 

3.1. 3:2 Resonance 

3.1.1. Dependence on planet m,ass, Mpi 

The first runs were undertaken with = 
2.5M0 and a = 60 AU. The nominal location of 
the planet for planetesimals at 60 AU to be in its 
3:2 resonance is api = 45.8 AU. Trapping probabil- 
ities were determined for planet masses of 1, 3, 10, 
30, 100, and SOOM©, and the results are plotted in 
Figure la. For a given planet mass, capture is as- 
sured as long as the migration rate is low enough. 
The capture probability is a strong function of the 
planet's migration rate, and the range of migration 
rates for which the capture probability is 0.1-0.9 
is relatively small. Such observations arc consis- 
tent with analytical predictions; e.g., a sharp de- 
pendence on the migration rate is predicted by 
autorcsonance theory (Friedland 2001), and cer- 
tain capture is expected in the adiabatic limit (i.e., 
when the migration rate is slow, see Appendix 
A) provided that the planetesimals' eccentricities 
are low enough (e.g., Henrard 1982; Henrard & 
Lemaitre 1983). The higher capture probabilities 
found with more massive planets for any given mi- 
gration rate are also expected, since their resonant 
gravitational perturbations are much stronger. 

A parametric fit to these trapping probabilities 



was performed that has the form: 



P = 



1 



Op/ 

ao.5 



n -1 



(3) 



where do. 5 is the migration rate for which half 
of the planetesimals are captured, and the pa- 
rameter 7 determines how fast the turnover is 
from a capture probability of 0.9 to 0.1 (e.g., 
(do.i - ao.9)/ao.5 = 9^/'^ - 'd'^/^). These fits are 
also plotted in Figure la. 

The parameters derived from these fits, do. 5 and 
7, are plotted in Figures lb and Ic to show how 
they vary with the mass of the migrating planet. 
A linear regression fit to the logarithm of these 
parameters is also shown on these plots. Thus it is 
found that the capture probability for the systems 
described in these runs can be approximated by 
equation (3) with the following parameters: 



= 0.153Mpi;3^<5±ooo4 



7 = 3.8M0-3«±0-0^ 



(4) 
(5) 



Equation (4) agrees well that the prediction from 
autorcsonance theory that the critical migration 
rate should scale with M^/^ (Friedland 2001). 
Equation (5) shows that lower mass planets have 
a larger range of migration rates resulting in in- 
termediate (0.1-0.9) capture probabihties. 

3.1.2. Dependence on semimajor axis, a 

Next the dependence of trapping probabilities 
on the planetesimals' semimajor axis was tested 
by performing a set of runs with M^, = 2.5Mq 
and Mpi = lOM®, and placing the planetesimals 
at semimajor axes of 30, (60), 100, and 150 AU. 
The results are plotted in Figure 2, where I have 
also plotted fits to the capture probabilities and 
to the derived parameters in the same manner as 
in §3.1.1. In this way the trapping probability is 
found to depend on a in the sense that do. 5 oc 
^-o.52±o.02_ dependence of 7 with a was found 
(7 a a0-00±0-05). 

3.1.3. Dependence on stellar mass, M^, 

Finally, the dependence of trapping probabil- 
ities on the stellar mass was tested. This was 
achieved with a set of runs with Mpi = lOM© 
and a = 60 AU, and trying different stellar masses 
of 0.5, 1.0, 1.5, (2.5), and 5.OM0. The results are 



4 



plotted in Figure 3, where fits to the capture prob- 
abilities and to the derived parameters are also 
plotted (see e.g., §§3.1.1 and 3.1.2). In this way 
the trapping probability is found to depend on Adi, 
in the sense that do.5 oc M-°-ss=^°-°^ A weak de- 
pendence of 7 oc M~^-^'^^-^^ was also found. 

3.1. 4-. Summary 

From the previous sections, it is evident that 
the trapping probability depends only on the di- 
mensionless parameters: 



e 



(6) 
(7) 



where fi determines the gravitational strength of 
the planet's resonance, and 9 is the ratio of the 
planet's migration rate to the planetesimal's or- 
bital velocity, which determines the angle at which 
the resonance is encountered. ^ Thus the proba- 
bility that planetosimals orbiting a star of mass 
Mi, (in Mq) at a somimajor axis a. (in AU) get 
captured into the 3:2 resonance of a planet of mass 
Mpi (in Mq) that is migrating at a rate dpi (in AU 
Myr~^) at a semimajor axis of (2/3)^/^a can be 
approximated by: 



P = 



1 + (X/i-"^) 



(8) 



Armed with this knowledge, the results of all 
the runs in §§3.1.1-3.1.3 were reanalysed to deter- 
mine A3:2, ^3:2, '^3:2, ^nd W3:2- A Icast squarcs 
fit to equation (8) was used to obtain best fit val- 
ues of: X3.,2 = 0.37, ^3:2 = 5.4, U3:2 = 1.37, and 
^^3:2 = 0.38. The errors in these parameters are es- 
timated to be ±0.02, 2.0, 0.02 and 0.1 respectively. 
This model for P was found to be accurate to 
about ±0.04 over the range of parameters tried 
in these runs. Translating back to the lexicon of 
equations (3)- (5), which can be compared with the 
parameters derived from Figs. (1-3): 



ao.5 
7 



2.7(Mp,/M0'-^V^*/«, 
5.4(Mp(/M0°-^^ 



(9) 
(10) 



Thus the results derived in §3.1.1 have not changed 
by the inclusion of the runs in §§3.1.2 and 3.1.3. 



^To allow you put the angle 9 in perspective, a migration 
rate of lAU Myr~^ for a IMq star corresponds to meeting 
the orbital velocity at 100 AU at an incident angle of 0^'3. 



3.2. 2:1 Resonance 

The same process was then repeated for the 
2:1 resonance, except that now that the scaling 
law is known, only the equivalent runs of §3.1.1 
needed to be performed. The results of runs with 
= 2.5Mq and a = 30 AU for different mass 
planets are shown in Figure 4a. ^ Fits of the 
form given in equation (3) were performed for each 
planet mass and fits to the variation of the derived 
parameters with Mpi were also performed. The 
resulting parameters were used to make initial es- 
timates of the equivalent parameters X2:i, ^2:1, 
U2:i, and V2:i in equation (8). A least squares fit 
to the capture probabilities then found these pa- 
rameters to be X2:l = 5.8, Y2:l = 4.3, U2:l = 1.40, 

and ?;2:i = 0.27, with respective errors estimated 
to be ±0.2, 2.0, 0.01, and 0.1. This model for P is 
shown on Figure 4a with dotted lines and is accu- 
rate to about ±0.025. 

The results for different resonances are summa- 
rized in Table 1. It is immediately obvious that 
the functional form of the capture probabilities for 
the two resonances are very similar and differ sig- 
nificantly only in the parameter X, which deter- 
mines the strength of the resonance. This numer- 
ical study finds that the 2:1 resonance is about 16 
times weaker than the 3:2 resonance. This is close 
to the factor of 12 predicted by autoresonance the- 
ory (Priedland 2001). 

3.3. 4:3 and 5:3 Resonances 

The same analysis as for §3.2 was then re- 
peated for the 4:3 and 5:3 resonances using runs 
with planetesimals orbiting initially 30 AU from a 
2.5M0 star. The results are plotted in Figures 4b 
and 4c, and the parameters derived from these re- 
sults given in Table 1. Again I find that for the 4:3 
resonance, the values of Y, u, and v are similar to 
those of the other first order resonances (i.e., those 
with q = 1). Also, X^:^ is close to that anticipated 
from autoresonance theory (i.e., that the 4:3 res- 
onance is ^ 1.6 times stronger than the 3:2 reso- 
nance). The values of Y, u and v for the second 



^While it might appear that integration times could be re- 
duced by achieving lower values of 9 by varying Mi, /a 
rather than reducing dpi, this is not the case, since the 

integrator chooses a stepsize that is oc ^Ji^^/Mi, and in- 
tegration length is oc Aa/dpi, where Ao oc a and dpi oc 
9sfMZl^. 
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order (q — 2) 5:3 resonance are somewhat differ- 
ent to those of the first order resonances; there is 
a steeper dependence of the strength of the reso- 
nance on the planet mass, and the range of migra- 
tion rates for intermediate trapping probabihties 
is greater for a given planet mass. In general sec- 
ond order resonances arc expected to be weaker 
than first order resonances, which is borne out by 
the higher value of X. 

3.4. Trapping Probability Summary 

The trapping probabilities derived in the last 
subsections are summarized in Figure 5. This 
shows, for the four resonances that were studied, 
the lines of 10, 50, and 90% trapping probability 
in terms of the migration parameters 6 and /U. For 
planctcsimals at a given scmimajor axis any given 
migration is defined by a single point on this Fig- 
ure. Thus this Figure quickly summarizes which 
resonances these planetesimals will end up in (if 
any). 

4. Distribution of Resonant Planetesimals 

Planetesimals that are trapped in resonance are 
not evenly distributed around star, rather they 
congregate at specific longitudes relative to the 
perturbing planet. This can be understood by 
first considering the geometry of a generic reso- 
nant orbit (§4.1), and then considering the action 
of resonant forces (§4.2). I use these in §4.3 to 
build up a model for resonant structure caused 
by planet migration, the parameters of which are 
determined from analysis of the numerical integra- 
tions described in §3. 

4.1. Resonant Geometry 

Each of the panels in Figure 6 shows the path 
that a planetesimal orbiting at a semimajor axis 
corresponding to a planet's p + q : p resonance 
(eq. [2]) takes when plotted in the frame co- 
rotating with the mean motion of the planet. 
While these orbits are elliptical in the inertial 
frame, there are two obvious features in the ro- 
tating frame: (i) when the planetesimal reaches 
pericenter it must be at one of p specific longi- 
tudes relative to the planet; (ii) the planetesimal 
spends more time at relative longitudes close to 
those at which it is at pericenter, than to those at 



which it is at apocenter. The former occurs be- 
cause by definition, the planetesimal's pericenter 
passages occur at longitudes relative to the planet 
incremented by q/p x 360° from the previous pas- 
sage, a pattern which repeats itself after p of the 
planetesimal's orbits. The latter occurs because 
when the planetesimal is close to its pericenter, 
the rate of change of its longitude is more closely 
matched to that of the planet. ^ 

While the pattern shown in the panels in Figure 
6 is unique, its orientation relative to the planet is 
not, since this is determined by the starting longi- 
tudes of the planet and planetesimal and the peri- 
center direction. This orientation can be defined 
using just one variable, the planetesimal's reso- 
nant argument: 

(/) = {p + q)Xj. - pXpi - qujr- (11) 

By definition, (j) remains constant while a planetes- 
imal is in resonance, since the increase due to 
is offset by the decrease due to Apj. The resonant 
argument has a specific physical meaning which 
will be described in §4.2, but from a geometrical 
point of view it can be rewritten thus: 

(j) = p[u}r - \pl{tperi)], (12) 

where Xpi{tperi) is the longitude of the planet when 
the planetesimal is at pericenter. In other words, 
(p/p defines the relative longitude when the plan- 
etesimal is at pericenter and so the orientation of 
the pattern. 

4.2. Resonant Forces 

The influence of a planet's gravity is to perturb 

the orbits of planetesimals in the system. These 
perturbations can be written down as a sum of 
many terms described by the planetesimal's dis- 
turbing function, R (Murray & Dermott 1999). 
Analysis of the disturbing function shows that 
these perturbations are particularly strong at the 
geometrical resonance locations (eq. [2]), where 
terms in which A,, and Ap; are combined in the 
form (p + q)Xr — pXpi are amplified. 

'In fact, even though the planetesimal orbits at a larger 
scmimajor axis, its longitude can change faster than that of 
the planet if the planetesimal's eccentricity is high enough. 
This results in backward motion in the rotating frame and 
causes the loops seen in Figure 6 for high eccentricity plan- 
etesimal orbits. 
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The resonances I will be discussing all involve 
terms including the angle (f) defined in equation 
(11). The relevant terms in the disturbing func- 
tion, assuming the orbits of the planet and plan- 
etesimal are coplanar and that the planet has a cir- 
cular orbit (i.e., the circular restricted three body 
problem), are: 



R = 



COS( 



a 

(13) 

where a = Upi/a, and i, /^j, and /j are coeffi- 
cients corresponding to the secular, resonant di- 
rect, and resonant indirect parts of the disturbing 
function respectively. The effect of these pertur- 
bations on the orbital elements of the planctesi- 
mal can be determined using Lagrange's planetary 
equations (e.g., Murray & Dermott 1999). In par- 
ticular, the variation in the planetesimal's semi- 
major axis and eccentricity are: 



_2(, + „(|5l),/^JM.e. 



a 



Mpi 



,9-1 



X (/d(a) + ^ ismc 



(14) 



(15) 



A simple analysis of the consequence of reso- 
nant forces using the circular restricted three body 
problem shows that they make the resonant argu- 
ment of a planetesimal librate about a fixed value 
(i.e., a sinusoidal oscillation; Murray & Dermott 
1999): 

= 0in + A0sin27r</t<^. (16) 

Indeed, a planetesimal is defined to be in reso- 
nance if its resonant argument is librating rather 
than circulating (i.e., a monotonic increase or de- 
crease). It is the fact that the resonant arguments 
of all planetesimals in the same resonance librate 
about the same value that causes their azimuthal 
distribution to be asymmetric. 

4.. 2.1. Libration Center without Migration 

The angle about which cj) librates, (j)m, can be 
understood by considering the geometry of res- 
onance (Peale 1976; Murray & Dermott 1999). 



Most of the perturbations to a planetesimal's orbit 
occur at conjunction (i.e., when the planet and the 
planetesimal are at the same longitude). Conjunc- 
tions which occur either at pericenter or apocenter 
confer no angular momentum to the planetesimal. 
However, conjunctions that occur before/after the 
planetesimal's apocenter (and after/before the 
pericenter) lead to a net decrease/increase in the 
planetesimal's angular momentum. This means 
that resonant forces push the conjunction towards 
apocenter. Since the resonant angle can also be 
written thus: 



= q{K - t^r), 



(17) 



where Ac is the longitude at which the planet and 
the planetesimal have a conjunction, apocentric 
libration is one at which 4>m/q = 180°. 

The argument outlined above is not valid how- 
ever for all resonances, just those with q = 1 
and p ^ 1. For a start, while conjimctions al- 
ways occur at the same location in the orbit of 
the planetesimal for a g = 1 resonance, conjunc- 
tions for q = 2 resonances occur at two loca- 
tions 180° apart around the planetesimal's orbit. 
Apocentric conjunctions would thus have to be fol- 
lowed by pericentric conjunctions. Since the forces 
of the pericentric conjunction would be stronger 
than those at apocenter, such libration is not sta- 
ble. Rather conjunctions for stable libration occur 
midway between pericenter and apocenter whence 
the forces from alternate conjunctions cancel, i.e., 
0^/2 = 90°. 

Also, while the resonant argument for the 2:1 
resonance librates about 180° for low eccentricity 
(e < 0.04) orbits, the libration center can take one 
of two values for higher eccentricity orbits, one 
with (f>m > 180° the other with 0^ < 180° (so- 
called asymmetric libration; Beauge 1994; Malho- 
tra 1996; Chiang & Jordan 2002). The physical 
explanation for this behavior is that perturbations 
to the orbits also occur when the planetesimal is 
at its pericenter, as well as when it is at conjunc- 
tion. Consider a 2:1 planetesimal which has a 
conjunction with the planet just before/after its 
apocenter. As already mentioned, the forces from 
the conjunction act to decrease/increase the plan- 
etesimal's angular momentum. Now consider the 
forces acting on the planetesimal as it continues 
along the rest of its orbit. At first these forces 
act to increase its angular momentum, since the 
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planet is in front of the planctcsimal, but when 
the planet's longitude is 180 — 360° in front of the 
planetesimal these forces decrease the planetesi- 
mal's angular momentum. Since the planetesi- 
mal spends more time at longitudes relative to the 
planet close to its pericenter (see Figure 6), these 
forces do not exactly cancel around the orbit and 
there is a net increase/decrease in the planetes- 
imal's angular momentum (for conjunctions be- 
fore/after apocenter). These forces are important 
for the p = 1 resonances for which every pericen- 
ter passage occurs at the same longitude relative 
to the planet. 

Stable libration for the 2:1 resonance occurs at 
a value of <!) which is < 180°/ > 180° for which 
the angular momentum loss/gain at conjunction is 
balanced by its gain/loss at pericenter. Including 
terms up to O(e^) in the expansion of the disturb- 
ing function for the circular restricted three body 
problem, the appropriate semimajor axis evolution 
at a Ri 2~^/^ is given by: 



a = — 




[3.38e sincj) 



+14.38e^sin2^- 2.52esin^]. (18) 

In this expression, the first two terms arc the di- 
rect terms of the p = q = I and p = q = 2 reso- 
nances respectively, and the third term is the in- 
direct term of the p = q = I resonance. The phys- 
ical interpretation of these terms is that the first 
two are from perturbations caused at conjunction, 
while the third term is from perturbations at peri- 
center. Setting d = gives the result that 

COS./)™ = -0.0298/e. (19) 

A similar result was found by Beauge (1994) by 
considering the phase space of the 2:1 resonance 
with a Hamiltonian model including harmonics up 

to second order. 

4-2.2. Libration Center with Migration 

The discussion of the stable libration center in 
§4.2.1 was based on the assumption that the res- 
onant forces confer no angular momentum to the 
planetesimal. This is not the case when the planet 
is migrating, since the planetesimal must be mi- 
grating out with the planet to maintain the reso- 
nance: 

/ _|_ \ 2/3 
amig = Clpl ( ) , (20) 



and so its angular momentum must be increas- 
ing. Given the discussion in §4.2.1, stable libra- 
tion should thus be offset to slightly higher values 
of (j) in the presence of migration. Setting equation 
(20) equal to the average rate of change of a due 
to resonant forces in the circular restricted three 
body problem (i.e., eq. [14] with (j) replaced with 
(pm) implies that 

sm(l>m(x-{e/n)e-i. (21) 

4-2.3. Eccentricity Evolution 

Another consequence of the planet's perturbing 
forces is to link the evolution of a planetcsimal's 
eccentricity to that of its semimajor axis. This 
means that the same resonant forces which cause 
a planetesimal's semimajor axis to increase while 
the planet's orbit is migrating out, also cause the 
planetesimal's eccentricity to increase. The eccen- 
tricity of a planetesimal which has migrated from 
an orbit with an eccentricity eg at a semimajor 
axis ao to one at a semimajor axis a due trapping 
in a p + q •- p resonance can be estimated from 
the circular restricted three body problem. Using 
equations (14) and (15), it can be shown that the 
change in a and e due to resonant forces are cor- 
related by the relationship da/de = 2(1 +p/q)ae- 
Thus the planetesimal's eccentricity can be esti- 
mated to be: 



e = Jel + ^^\na/ao. (22) 

V p + q 

4.3. Numerical Results 

Based on the discussion of §§4.1 and 4.2, it is 
easy to see that the spatial distribution of a popu- 
lation of planetesimals that have been trapped in 
resonance by a migrating planet can be defined by 
the distribution of their orbital parameters (see 
also §5.1). Orbital distributions sufficient to do 
this are derived in this section for the planetesi- 
mals in the migration simulations presented in §3. 

First the output of these simulations were used 
to determine the distributions of the longitudes 
of the resonant planetesimals as well as the evo- 
lution of their eccentricities while in resonance. 
Then additional runs were performed to determine 
the libration parameters of the resonant planetes- 
imals. These used the output of the original simu- 
lations (i.e., the instantaneous orbital parameters 
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of the planctcsimals and planet) as their starting 
point and continued the evolution for just a few li- 
bration periods (see Appendix A) sampled at 500 
timcstcps. The libration parameters 0„j and 
for each planetesimal in the run that was trapped 
in resonance were then measured by a fit to the 
evolution of its appropriate resonant argument us- 
ing equation (16). 

In order to determine the mean parameters for 
the ensemble of planctcsimals that were trapped 
in resonance, a histogram of these parameters was 
displayed. The mean libration center of these 
planetesimals, (^m), was determined, and a gaus- 
sian fit performed to the distribution of libration 
amplitudes. This resulted in the best fit mean li- 
bration amplitude, (Acj)) and standard deviation 
of the distribution of libration amplitudes, cta.^; 
a gaussian distribution always provided a decent 
approximation for the distribution of libration am- 
plitudes. 

4.3.1. Distribution of Longitudes 

It had originally been anticipated that the dis- 
tribution of the longitudes of resonant planetesi- 
mals in the runs would be random, since they were 
started with random longitudes. This was not the 
case J a typical (but not unique) feature of the lon- 
gitude distributions obtained was a double-peak 
superimposed on a random distribution. This was 
found to be an artifact of the initial conditions, 
caused by all planetesimals being started at ex- 
actly the same scmimajor axis. When an addi- 
tional population of planetesimals at a different 
semimajor axis were introduced into a run, these 
were also captured into resonance and their lon- 
gitude distribution was also double-peaked. How- 
ever, at any one time, the peaks for the popula- 
tions which originated at different scmimajor axes 
occurred at different longitudes. Further runs with 
planetesimals started with semimajor axes ran- 
domly chosen within a range, and their longitude 
distribution was indeed found to be random. Since 
the planetesimals that end up in resonance origi- 
nate from a range of semimajor axes, their longi- 
tudes can be assumed to be random after resonant 
trapping. This means that the instantaneous dis- 
tribution of a population of trapped planetesimals 
which all have the same resonant argument and ec- 
centricity would look like that shown in Figure 6, 
since this figure shows paths of such planetesimals 



at regular increments in their longitudes. 

4.3.2. Eccentricity Increase 

The eccentricity increase for planetesimals that 

are trapped in resonance was found to be well ap- 
proximated by equation (22) for all of the runs. 

4.3.3. Libration Centers 
3:2, 4:3 and 5:3 Resonances 

As expected on physical groimds (§4.2.1), for the 
three resonances 3:2, 4:3, and 5:3, the libration 
centers were found to tend towards 180° for low 
planet migration rates. The way {(pm) responded 
to increasing the planet's migration rate (§4.2.2) 
and varying the other parameters was tested using 
the sets of runs for the 3:2 resonance which varied 
the planet mass, planetesimal semimajor axis and 
stellar mass. These runs showed that {^m) only 
depends on 6/^. 

Figure 7 shows the deviation of {(j>m) from 
180° for each of the resonances as a function of 
6/i.j,. These show an approximately linear cor- 
relation, as expected from the circular restricted 
three body problem (eq. [21]), but with a turnover 
for high However, only a very small vari- 

ation of {(pm) was found in the course of migra- 
tion as e increases, serving as a caution against 
using simple analytical solutions of the circular 
restricted three body problem. Fits for each of 
the resonances were performed having the form 
{(l)m) - 180° = A{e/iJ,) + B{e/n)'^. The results are: 

(</.„5,3) - 180° = 5.8(^//x), (23) 
((^™3,2) - 180° = 7.5(0/m) - O.23(0//z)2(24) 
(<^„4,3) - 180° = 5.3{9/tx)- 0.12(6/ fi)\25) 

and these fits are also shown on Figure 7. 

2:l(u) and 2:1(1) Resonances 

Again as expected on physical grounds (§4.2.1), 
the libration center of a planetesimal in the 2:1 
resonance was found to change during the course 
of the migration as its eccentricity increased. Also, 
the evolution of its libration center was found to 
take one of two courses: one in which increased 
during the migration such that it was always > 
180°, and the other in which (f)^ decreased and so 
was always < 180°. From now on I will refer to 
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these as two separate resonances, the 2:l(u) and 
2:1(1) resonances, respectively. 

The mean libration centers of the resonant 
planetesimals in the runs shown in Figure 4a are 
plotted in Figure 8; planetesimals were separated 
into the 2:l(u) and 2:1(1) resonances and their li- 
bration centers measured at four epochs through- 
out the migration. The libration centers in these 
runs were found to depend almost entirely on the 
planctesimal's cicccntricity, with any variation due 
to fj, 01 9 barely discernible. Neither was any sig- 
nificant difference found between the 2:l(u) and 
2:1(1) resonances, and a simple parameterised fit 
was performed of the form cos {(l>m) = A + B/e. 
The result is: 

cos {(I)m2..i) = 0.39-0.061/e, (26) 

and these lines are shown on Figure 8 with the 
additional constraint that cos((^m2 i) ^ ~1- The 
line derived from the circular restricted three body 
problem (eq. [19]) is also shown on this figure. 

4-3. 4- Trapping Probability for 2:1(1) Resonance 

The trapping probabilities derived in §3.2 do 
not take into account whether the planetesimal 
ends up in the 2:l(u) or 2:1(1) resonance, since 
the semimajor axis and eccentricity evolution is 
the same for both resonances. However, it was 
recently reported that trapping probabilities are 
not the same for the two resonances (Chiang & 
Jordan 2002). Thus I reanalysed the trapping 
probability runs for the 2:1 resonance to determine 
which of the resonances was being populated; the 
results are shown in Figure 9, where I have plot- 
ted P2:i(;): the trapping probability for the 2:1(1) 
resonance. It was found that there are always less 
planetesimals trapped in the 2:1(1) resonance than 
in the 2:l(u), except in the limit where the planet's 
migration rate tends to zero at which point the 
two resonances axe equally populated. ^ As the 
migration rate is increased, while trapping into 
the 2:1 resonance is still 100%, -P2:i(/) decreases 
oc ^0.5^-0.25 ^^]\ reaches zero. This behavior 

^ Since the original runs were designed to have trapping prob- 
abilities close to 50%, additional runs had to be carried out 
to ascertain how -P2:i(!) varies as planet migration rate is 
reduced when trapping is 100% for the 2:1 resonance as a 
whole. 



can be approximated by: 

Pi-.m = 0.5 - 0.856i°- (27) 

Unusually, as the migration rate is increased 

further to the point where trapping into the 2:1 
resonance begins to decrease, trapping probabil- 
ities for the 2:1(1) resonance increase again such 
that a few per cent can become trapped. For mod- 
eling purposes this behavior was approximated us- 
ing the following functions: If > 0.09 then 
the capture probabilities given in equation (27) 
should be increased by an amount 

dP2:i(i) = 0.11 - 0.486'/x-^■^^ (28) 

assuming that dP2:i(i) > 0. The complete fits are 
shown in Figure 9, and the line delineating mi- 
grations for which P2:i{i) ~ 1/3 (i.e.. for which the 
2:l(u) resonance is twice as populated as the 2:1(1) 
resonance) is shown in Figure 5 (i.e., 6 w 0.038y^ 
from eq. [27]). 

The reason for this asymmetry remains unclear, 
however a clue to its origin comes from the physi- 
cal interpretation of the origin of the resonance. 
As the planctesimal's eccentricity increases, its 
resonant argument initially librates about 180°. 
Once its eccentricity reaches ~ 0.03 the libration 
center must either increase or decrease to balance 
the angular momentum imparted to it at conjunc- 
tion and pericenter (§4.2.1). Nothing in the argu- 
ment presented so far has hinted at either of the 
resonances being stronger than the other. How- 
ever, because of the migration, the stable libra- 
tion is not exactly at (pm = 180°, but at a slightly 
higher value (§4.2.2), even if this offset is too small 
to detect in these runs (§4.3.3). The fact that (p is 
more often (if not always) > 180°, may make the 
2:l(u) resonance the more likely path. This would 
tie in with a higher asymmetry expected for higher 
migration rates for which the offset of the libration 
center from 180° is higher (§4.2.2). 

4-3.5. Libration Amplitude Distributions 

The distribution of libration amplitudes was al- 
ways found to be fairly broad with <ta<i> in the 
range 5 — 20° . No significant correlation could be 
found of aA4, with either /j. or 6 for any of the res- 
onances. Thus for modelling purposes these dis- 
tributions were assumed to have values of 

(TA4>2a = 10°, (29) 
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CA03,2 = 15°, 
CrA04:3 = 15°, 



(30) 
(31) 
(32) 



which are the means of all the runs performed for 
each resonance. 

However, the mean libration amplitudes, (A^), 
which arc plotted in Figure 10, were found to 
vary in a systematic manner determined by the 
dimensionless parameters fx and 6. In particular 
the mean libration amplitudes depend only on the 
combination O/fj,^, so that libration amplitudes 
are higher for faster migrations and more mas- 
sive planets. For the 2:1 resonance the libration 
amplitude was also found to decrease during the 
migration as the planetesimals' eccentricities are 
increased. No significant variation of (A0) dur- 
ing the migration was found for the 5:3, 3:2, and 
4:3 resonances, and only {A(f)) at the end of their 
runs were considered. Also, no significant differ- 
ence was found between the 2:l(u) and 2:1(1) res- 
onances, and so their results were considered to- 
gether as showing the variation for the 2:1 reso- 
nance. 

Least squares fits were performed to these re- 
sults of the form (Acj)) = A + B{9/ fi^)/e'^ , where 
D = except for the 2:1 resonance. The results 
of these fits are: 

{Act>2:i) = l.l + 42.2(0//xi-24)/eO-42, (33) 

(A(/)5:3) = 13.2 + 290.0(e//z^-3''), (34) 

(A</,3:2) = 9.2+11.2(^//xi-27), (35) 

(A04:3) = 5.0 + 7.9(^/^1-27), (36) 

and these fits are also shown in Figure 10. 

4-3.6. Other Parameters 

Other parameters were also derived from these 
runs which arc interesting from a celestial me- 
chanics point of view, but which are not directly 
relevant for the structure of a planetesimal disk. 
These are described in Appendix A. 

5. Debris Disk Structure 

The results given in §§3 and 4 can be used to 
predict the spatial distribution of planetesimals at 
the end of any given migration defined by fj, and 9 
(eqs. [6] and [7]), and in §5.1 a numerical scheme 
is described which does just that. Perhaps more 



importantly, the converse is also true: an observed 
spatial distribution of planetesimals can be used to 
constrain the migration which caused that distri- 
bution. To help with such an interpretation, some 
generalizations about the kind of structures that 
result from different migrations are given in §5.2. 
In §5.3 this model is applied to observations of 
Vega's debris disk which shows that it is possible 
to set quite tight constraints on the migration his- 
tory of this system. §5.4 discusses the limitations 
of the model and future developments. 

5.1. Numerical Model 

The model is completely defined by the follow- 
ing input parameters: 

• Planetesimals: their initial distribution is 



defined by a,, 



and (5, where the 



number of planetesimals in the semimajor 
axis range a to a + da is oc a^da; in the min- 
imum mass solar nebula model 6 = —0.5. 

• Planet: has a mass Mpi and a migration 
defined by Upi., Upi^, and dpi, which are its 
initial and final semimajor axes and its mi- 
gration rate, respectively. 

• Star: has a mass M*. 

A number of planetesimals, Npp, arc then dis- 
tributed in semimajor axis randomly according 
to the prescription above, with eccentricities also 
randomly chosen between and 0.01; Npp has to 
be sufficient to define the distribution of eccentric- 
ities of resonant planetesimals in the final distri- 
bution, and is normally set at 400. For each plan- 
etesimal, the passing of each of the planet's reso- 
nances is considered in the order they encounter 
the planetesimal. If a random number in the range 
0-1 is less than the trapping probability defined 
by equation (8) and Table 1, then the planetes- 
imal is assumed to become trapped in that res- 
onance; for planetesimals that are trapped in the 
2:1 resonance, the probability that this is the 2:1(1) 
resonance is also determined from equations (27) 
and (28). Naturally no more of the resonances are 
then considered, and this planetesimal's eccentric- 
ity and semimajor axis are assumed to increase 
according to equations (2) and (22). 

Each of the resonant planetesimals is assumed 
to be representative of 9600(p + g) more. For each 
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planctcsimal, 400 values of (j) were chosen by first 
using equations (23) - (26) to determine tlie appro- 
priate libration center. Then 20 libration ampli- 
tudes were chosen at random from the appropriate 
gaussian distribution defined by equations (29) - 
(36), and for each libration amplitude 20 values of 
(j) were chosen from equation (16) with values of 
t/t(j, evenly distributed between and 1. For each 
of these values of (f), 24(p + q) planetesimals were 
assigned evenly spaced longitudes so that their 
spatial distribution matched the patterns shown in 
Figure 6 (with an appropriate orientation defined 
by (j)). The resulting number density distributions 
for each of the resonant planetesimals was then 
normalised to unity (by dividing by 9600 (p + q))- 
The semimajor axes and eccentricities of non- 
resonant planetesimals remain unchanged during 
the migration unless they reach the planet's reso- 
nance overlap region: the region 

\a/api - 1| < 1.3(3 x lO'V)^^^ (37) 

is chaotic and planetesimals entering this region 

would be scattered out of the system on short 
timescales (Wisdom 1980); such planetesimals are 
removed from the model. Each of the remain- 
ing non-resonant planetesimals was assumed to be 
representative of 400 more with even distributions 
of longitude and longitude of pericenter, and their 
resulting number density distributions normalised 
to unity (by dividing by 400). 

5.2. Generic Structures 

First consider the distribution of planetesimals 
ignoring both the change in libration center due 
to migration and the amplitude of the libration, 
i.e. such that (j){t) = (j)mifi = 0). The spatial 
distribution of such planetesimals would look like 
that shown in Figure 6. That is, planetesimals 
in the 3:2 resonance would be strongly concen- 
trated at ±90° longitude relative to the planet, 
while those in the 4:3 and 5:3 resonances would 
be concentrated ±60° and 180° from the planet, 
and those in the 2:l(u) resonance would be con- 
centrated 103 — 79° behind the planet (for eccen- 
tricities 0.1-0.3) and most of those in the 2:1(1) 
would be found 103 — 79° in front of the planet. 
The strength and physical size of the concentra- 
tions of a population of planetesimals that are in 
a given resonance depend on the distribution of 
their eccentricities. 



The fraction of planetesimals that end up in a 
particular resonance is determined not only by the 
H and 6 of the migration, but also by the initial 
distribution of its planctcsimal population, as well 
as the extent of the migration. To help visualize 
what the outcome of any given migration would 
be. Figure 11 and Table 2 summarize the dynam- 
ical structures resulting from migrations that are 
located in the seven different zones in the ijl—6 plot 
of Figure 5. The boundaries between the zones 
are taken as the lines of 50% trapping probabil- 
ity for the different resonances, as well as the line 
for which twice as many planetesimals arc in the 
2:l(u) resonances compared with the 2:1(1) reso- 
nance. Clearly these boundaries are not so dis- 
tinct, although the areas covered by 10-90% trap- 
ping probabilities are relatively small on this plot. 
The application of these zones will become clearer 
in §5.3. 

Consider now the offset in libration center due 
to migration. This does not affect the spatial 
distribution of 2:1 planetesimals, but results in a 
rotation of the pattern for the other resonances 
shown in Figure 6 by an amount (</>„ — 180°)/p. 
This is plotted in Figure 12 for migrations which 
result in trapping probabilities of 10, 50, and 90%. 
Since the turnover is not well characterized for 
high values of // the line is not extrapolated above 
the highest value of /z in the runs. The rotation is 
always negligible for the 5:3 resonance, and is usu- 
ally small, < 30°, for the 3:2 and 4:3 resonances. 
In other words, this rotation would only be no- 
ticeable for high mass planets (/U ~ 100) migrat- 
ing close to the limit where trapping probabilities 
are < 50%. Further, this rotation is only valid 
while the planet is migrating. Simulations similar 
to those already described were performed with 
the planet migration turned off. In this instance 
the libration was about 180° (apart from the libra- 
tion centers of the 2:1 resonances which were un- 
changed). Since the rotation is small it is included 
in the model for consistency, but any model which 
relies on this rotation must consider the probabil- 
ity of our witnessing a system mid-migration. 

Consider now the effect of the libration of (p 
due to a non-zero A(/). Since the oscillation is si- 
nusoidal, the distribution of (j) is not peakc^d at 4>m, 
but actually has two peaks at (j)m ± A0. In princi- 
ple, if A(/) is big enough, the planetesimal distribu- 
tion will peak at 2p rather than p longitudes rel- 
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ativc to the planet (Chiang & Jordan 2003). The 
maximum rotation of the pattern from the orien- 
tation defined by 0m (i-e. that shown in Figure 6) 
is ztA(j)/p, and Figure 13 shows A(f)/p for all res- 
onances as a function of /z for the values of for 
which trapping probabilities are 10, 50, and 90%. 
For the 2:1 resonance, the libration amplitude de- 
creases throughout the migration, and so is plotted 
at two reference eccentricities, e = 0.03 and 0.3. 
Given the physical size of the clumps, these libra- 
tion amplitudes are relatively small, even when the 
migration is fast enough to result in low trapping 
probabilities. Thus libration in this model results 
only in a slight azimuthal smearing of each clump. 

5.3. Application to Vega 

Vega is a nearby (7.8 pc) main sequence AO 
star (M^ « 2.5Mq) with an age of ~ 350 Myr 
(Song et al. 2001). Its emission spectrum exhibits 
a strong excess above the level of the photosphere 
at wavelengths longward of 12 /xm (Aumann et al. 
1984). This excess originates in dust grains or- 
biting the star that are continually replenished by 
the destruction of larger planetesimals in its debris 
disk. Imaging at submillimeter wavelengths shows 
the morphology of the excess emission (Holland et 
al. 1998; see Fig. 14a) down to a resolution of 
14" (FWHM): the lowest contours are nearly cir- 
cularly symmetric extending to ^ 24", indicating 
that the disk is being observed pole-on, but the 
disk's structure is dominated by an emission peak 
^ 9" to the northeast of the star; the highest con- 
tours are also elongated in the southwest direction. 

The structure of the disk was recently mapped 
with even higher spatial resolution (^ 3") using 
millimeter interferometry (Koerner, Sargent & Os- 
troff 2001; Wilner et al. 2002). While these obser- 
vations were not sensitive to the larger scale struc- 
ture, they did show that a significant fraction of 
the millimeter emission could be resolved into two 
clumps, one 9.5" northeast of the star, the other 
8" southwest of the star, and that the northeast 
clump is brighter than that southwest. Thus it 
appears that the submillimeter images arc best in- 
terpreted as emission from an extended disk which 
is dominated by two roughly equidistant clumps of 
unequal brightness on opposite sides of the star. 

Based on the discussion of §5.2, the migration 
zone (Fig. 11) causing this structure can be nar- 
rowed down to zone Di: Two clumps of equal 



brightness could have been explained by planetes- 
imals trapped in the 3:2 resonance (zone C), or 
in equal numbers into the 2:l(u) and 2:1(1) reso- 
nances (zones Dii or Eii), but the asymmetry in- 
dicates that one of the clumps is overpopulated, 
pointing to migration zone Di or Ei. Further con- 
straints are set by the lack of evidence for an ad- 
ditional 3 clump pattern rotated relative to the 
two clump pattern. While 3 clumps from the 4:3 
resonance arc inevitable, zone Ei is ruled out as 
trapping of planetesimals into the 5:3 resonance 
in this zone means there are less available to be 
trapped into the 2:l(u) resonance. In theory the 
outcome of migrations anywhere within zone Di 
will be similar and are not constrained by this 
model, however the fuzzy edges of the zones mean 
that further constraints are possible. Here I set 
the mass of the planet to be the same as that of 
Neptune, Mpi = 17.2M0 = 6.9), which means 
that the migration rate must be close to 0.5 AU 
Myr~^ 

The remaining parameters were then con- 
strained by a best fit to the submillimeter disk 
image presented in Holland et al. (1998) and re- 
produced in Figure 14a: the only variables were 
o-max, dpi-, CLpij, and hpi, since 5 was fixed at -0.5 
and amin was set at ap; . . The observing geometry 
was assumed to be face-on and the spatial dis- 
tribution of planetesimals was converted into an 
image of the dust emission resulting from the de- 
struction of those planetesimals using the follow- 
ing assumptions: (i) that the spatial distribution 
of the dust exactly follows the distribution of the 
parent bodies (i.e., so that the cross-sectional area 
of emitting dust is proportional to the number of 
planetesimals); (ii) that the grains' submillimeter 
emission is oc which is a good approxima- 

tion for the large grains which contribute to a 
disk's submillimeter emission. 

The resulting best fit is shown in Figure 14b, 
and was performed by comparison of the con- 
tours of the two images; the orbital and spatial 
distributions of the underlying planctcsimal pop- 
ulation arc shown in Figure 15. The best fit 
parameters were found to be amax = 140 ± 15 
AU, api. = 40 ± 10 AU, Upij = 65 ± 5 AU, and 
hpi = 0.45 ± 0.1 AU Myr"'^ (corresponding to 
6* = 1.8 - 3.4 at 40-140 AU and a total migration 
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time of 56 Myr^). In the final distribution 5.1% 
of the planetesimals are trapped in the 4:3 reso- 
nance, 22.5% in the 3:2 resonance, 0.4% in the 5:3 
resonance, 18.6% in the 2:l(u) resonance, 0.7% in 
the 2:1(1) resonance, 41.0% remain in nonresonant 
orbits, while 11.6% are ejected by resonance over- 
lap. The outer edge of the disk, amax, was con- 
strained to within ~ 15 AU by a fit to the lowest 
contours in the images. The final planet location, 
ttpi^ , was constrained to within ^ 5 AU by the ra- 
dial offset of the northeast clump. The contrast of 
the clumps as well as their morphology were then 
used to constrain the initial location of the planet, 
Upi. , and the migration rate, dpi . Since the best fit 
has a migration rate for which trapping into the 
2:1 resonance is < 100%, and into the 5:3 reso- 
nance is > 0%, small changes in dpi of the order 
±0.02 AU Myr~^ resulted in large changes in the 
model structure. However, dpi could not be so well 
constrained, because the extent of the migration, 
and so Upi., also has a large effect on the model 
structure, since this determines both the fraction 
of planetesimals in each resonance as well as their 
eccentricity distributions. The errors given above 
are not formal errors, but approximate limits for 
which acceptable fits are possible to the structure 
of Figure 14a. 

This demonstrates that this model can fit the 
observations very well and in doing so constrains 
important parameters regarding the evolution of 
this system. The implications of this model are 
as follows. First of all the whole pattern is ex- 
pected to orbit the star with the same period as 
the planet. As this is predicted to be at 65 AU, 
the orbital period is 330 yr. Since the clumps 
are ~ 9" from the star, their motion would be 
0'.'17 yr^^. With an absolute pointing miccrtainty 
of ±2" from the James Clerk Maxwell Telescope, 
such absolute motion would not be detectable for 
several decades in submillimcter images with cur- 
rent technology. However, it would be easier to 
detect the relative motion of the two clumps, since 
the position angle of this structure should vary at 
1.1° yr~^. The direction of the orbital motion pre- 
sented in Figure 14 is anticlockwise. However this 

® Given the age of the system, this implies that the migra- 
tion is now finished, thus it is important to point out that 
the rotation of the resonant structure in the model due to 
migration (i.e., 4>rn — 180°: SS4.2.2 and 5.2) is very small, 
< 2° , and so its effect on the derived structure is negligible. 



solution is not unique, since the model is nearly 
symmetrical about the line joining the two clumps, 
and I could also have modeled the system with 
clockwise orbital motion, in which case the planet 
would be currently north of the star. Such a model 
is shown in Figure 16; there are only slight dif- 
ferences in the disk structure compared with Fig- 
ure 14b. The important point is that the brighter 
clump, corresponding to the 2:l(u) resonant plan- 
etesimals, trails some 75° degrees in longitude be- 
hind the planet. 

Another implication of the model is that the 
emission distribution is much more intricate than 
that detectable with a 14" beam. It is therefore 
interesting that in interferometric images of Vega 
(Koerner et al. 2001), the northeast clump splits 
into three, with two lower level clumps straddling 
a brighter clump; the brighter clump in the model 
presented in this paper is expected to be extended 
in longitude (see Figure 15b). This model also 
makes the prediction that a lower level three clump 
pattern exists from planetesimals in the 4:3 reso- 
nance. While two of the clumps almost blend with 
the two brighter clumps from the 3:2 and 2:1 plan- 
etesimals, a faint clump is expected on the oppo- 
site side of the star from the planet. If the presence 
of such a clump could be inferred from higher reso- 
lution and sensitivity observations of this disk, the 
location of the planet and its direction of motion 
could be unambiguously determined. 

The mass and migration rate of the planet have 
been constrained to lie at the lower edge of zone Di 
of Figure 1 1 , although an additional constraint is 
that the migration rate must be > 0.07 AU Myr^^ 
for the 25 AU migration to have been completed 
over the age of the system. Thus the planet mass 
and migration rate I chose is not unique and must 
be determined by a study of the origin of the mi- 
gration rate. However, it is noteworthy that in 
Hahn & Malhotra (1999), Neptune was found to 
migrate from 23-30 AU over 50 Myr in their model 
with an initial planetesimal disk mass of 5OM0, 
a migration rate not too dissimilar to that cho- 
sen for this model. In other words, the mass and 
migration rate of the planet causing the struc- 
ture of Vega's disk could be similar to the mass 
and migration rate of Neptune in our own sys- 
tem. This similarity to the solar system is con- 
trary to previous models for Vega's disk structure, 
which was interpreted as dust grains migrating in 
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toward the star due to Poynting-Robcrtson drag 
that get trapped into the 2:1 and 3:1 resonances 
of a 3 Jupiter mass planet on an orbit with an 
eccentricity of 0.6 (Wilner et al. 2002). 

As for the possibihty of directly detecting this 
planet, observational constraints have recently 
been set on the presence of planets around Vega, 
with an upper limit of ~ 33OOM0 at 10" from the 
star (Metchev, Hillenbrand & White 2003). Given 
its proximity to the bright star Vega, it would 
be very difficult to detect this planet with cur- 
rent technology if it is indeed as small as 17 M^. 
A larger planet, with a necessarily much faster 
migration rate, may be detectable. 

5.4. Caveats 

Despite the complexity of the modc;l, it is still 
just a first step toward a complete explanation of 
the structures which could be caused in extrasolar 
systems by planet migration. 

For a start, the trapping probabilities were de- 
rived for migration through a dynamically cold 
disk (e < 0.01). The planetesimals' eccentricities 
have a significant infiuence on trapping probabil- 
ities, which decrease as the eccentricities are in- 
creased, and also do not reach a maximum at 100% 
trapping probability for low migration rates if the 
eccentricity is above a certain threshold (Borderies 
& Goldreich 1984: Melita & Brunini 2000). Most 
importantly, the trapping probabilities of differ- 
ent resonances are affected in different ways, with 
higher order (q > 1) resonances becoming more 
populated relative to lower order resonances for 
higher eccentricities. Indeed Chiang et al. (2003) 
propose that the Kuiper belt must have been dy- 
namically hot (e w 0.2) when the migration of 
Neptune occurred to explain the presence of three 
KBOs in the 5:2 resonance given that only seven 
objects have been discovered in Neptune's 2:1 res- 
onance. It is difficult to predict an appropriate 
eccentricity distribution for the residual planetcs- 
imal disk. While planetesimals are thought to 
form in a dynamically cold disk, the sweeping 
of a planet's resonances and its secular perturba- 
tions, as well as scattering of planetesimals from 
closer to the star, can all contribute to increasing 
the average eccentricity of the planetesimal disk. 
Certainly future models will have to consider the 
possibility of high eccentricities, but in doing so 
will become more complicated, since trapping into 



higher order resonances will also have to be con- 
sidered. 

Another omission of the model is stochastic ef- 
fects. Migration in the models of Hahn & Malho- 
tra (1999) is not smooth, but shows large jumps. 
This is because the residual planetesimals used in 
their simulation had to bo large to give reasonable 
integration times. If the most massive residual 
planetesimals causing migration are smaller than 
a certain limit, then stochastic effects can be ig- 
nored (Chiang et al. 2003), otherwise they must be 
characterized (Zhou et al. 2002). One reason for 
anticipating that the residual planetesimals caus- 
ing the migration are much smaller than the planet 
would be if the planet formed much closer to the 
star than the planetesimals. This could be the 
case if the planet was flung out in a gravitational 
interaction with other giant planets (Thommes et 
al. 2002) or migrated out to a more distant loca- 
tion. 

Possibly the most worrying aspect of the model 
is the translation of the distribution of planetesi- 
mals to the distribution of dust. These were as- 
sumed to be identical, but this may not be ex- 
actly true. For a start the radial location of the 
resonances are different for small dust, since they 
orbit the star slower than larger grains due to ra- 
diation pressure (Wyatt et al. 1999). Also, these 
particles have different orbits from their parents, 
not only because of radiation pressure, but also 
due to the velocity given to them during the colli- 
sion. The effect of P-R drag and subsequent col- 
lisions should also be taken into account. These 
are important issues, but ones which merit a pa- 
per in themselves. However, in defense of this as- 
sumption I will point out two things. First, the 
dust which contributes to the submillimeter im- 
ages is expected to be large, since small grains 
emit very inefficiently at such long wavelengths 
(Wyatt & Dent 2002). Second, most resonant 
planetesimals are on planet-crossing orbits (i.e., 
ares(l — fires) < Opj), and the only reason such or- 
bits are stable is because resonant forces prevent a 
close approach (§4). Small dust originating from 
such planetesimals, but which is no longer in res- 
onance, would therefore be expected to be short- 
lived, since without the protection of the resonance 
a close approach to the planet is possible resulting 
in scattering out of the system. Such scattering is 
confirmed in more recent runs following the evolu- 
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tion of the dust particles' orbits (Wyatt, in prep.). 

Another concern is that the dust distribution 
may be affected by random events in which single 
massive planetesimals are disrupted. Such events 
may have been witnessed in the structure of the 
zodiacal cloud (Dermott et al. 2002; Nesvorny et 
al. 2003). However, Wyatt & Dent (2002) con- 
cluded that inidividual collisions are unlikely to 
be a significant source of structure in the Fomal- 
haut disk (which also contains a clump, Holland 
et al. 2003), since this disk is so massive that 
only a collision between two planetesimals at least 
1400 km in diameter could affect its structure at 
the level observed, and too few such planetesimals 
can coexist in the disk. Since Vega's disk, and its 
clumps, have a similar mass to those of Fomalhaut 
(Holland et al. 1998), I conclude that random col- 
lisions are also unlikely to be the cause of Vega's 
clumpy disk structure. 

While I prefer to leave the discussion of the 
origin of the planet's migration to a later paper, 
it is worth pointing out that only certain values 
of 9 will be possible for a given planet mass /i 
and initial planetesimal disk mass (Hahn & Mal- 
hotra 1999). In other words, the planet's mass 
and migration rate, as well as the initial mass of 
Vega's planetesimal disk, could be much better 
constrained with a model incorporating both mi- 
gration and resonant trapping. Also, in the models 
of the migration of the solar system's planets, Nep- 
tune's migration is outward only because of the 
existence of the massive planets interior to its or- 
bit. This is because the planetesimals scattered in- 
wards by Neptune (causing its outward migration) 
may not retiirn for a subsequent scattering event 
(causing inward motion) due to their interaction 
with the closer in planets, whereas the converse is 
not true. Thus an outward migration may be said 
to be indicative of at least two planets. However, 
Ida et al. (2000) presented an analytical model for 
planet migration claiming that once outward mi- 
gration has started, it is self-sustaining, since there 
are automatically less planetesimals interior to the 
planet's orbit. 

6. Conclusion 

I have presented a model describing the conse- 
quence of the outward migration of a planet on 
the dynamical and spatial structure of a planetes- 



imal disk residing outside the planet's orbit. In §3 
numerical simulations were used to show how trap- 
ping probabilities into the 4:3, 3:2, 5:3, and 2:1 res- 
onances can be estimated to within 5% given just 
two parameters = Mpi/M^, and 9 = apis/a/M^. 
Resonant forces cause a planetesimal's resonant 
argument, (p, to librate and in §4 physical argu- 
ments were used to explain what azimuthal struc- 
ture is expected in the distribution of planetesi- 
mals clS Qb result of this libration: planetesimals 
trapped in the 4:3 and 5:3 resonances are concen- 
trated at 60, 180, and 300° longitude relative to 
the planet, those in the 3:2 resonance are concen- 
trated at ±90° longitude, and those in the 2:1 res- 
onance are concentrated in two clumps associated 
with the 2:l(u) resonance with a concentration at 
relative longitude of ~ 285°, and the 2:1(1) reso- 
nance with a concentration at ~ 75°. Also in §4, 
numerical simulations were used to show how the 
angle about which </> librates, as well as the ampli- 
tude of that libration, are affected by the migra- 
tion parameters. These simulations also charac- 
terized the overpopulation of the 2:l(u) resonance 
relative to the 2:1(1) for different migrations. In §5 
the numerical results were used to derive a numer- 
ical scheme to predict the spatial distribution of 
planetesimals resulting from a migration defined 
by /i and 9. It was shown that the dynamical 
structure of a post-migration disk can have one 
of seven states depending on the /x and 9 of the 
migration. 

Application of the model to the structure 
of Vega's debris disk presented in Holland et 
al. (1998) shows that its two clumps of unequal 
brightness can be explained by the migration of 
a Neptune mass planet from 40-65 AU over ~ 56 
Myr through a planetesimal disk initially extend- 
ing from 40-140 AU. The two clumps are caused 
by planetesimals trapped in the 3:2 and 2:1 res- 
onances, and the brightness asymmetry is caused 
by an overabundance of planetesimals in the 2:l(u) 
resonance relative to the 2:1(1) resonance. While 
the extent of the planet's migration is well con- 
strained by the brightness of the clumps, its mass 
and migration rate are not unique, although they 
are constrained to lie within certain ranges defined 
by zone Di in Figure 11. Further constraints on 
the planet's mass and migration rate, as well as 
on the mass of the planetesimal disk, would be 
possible by modeling the origin of the planet's mi- 
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gration in angular momentum exchange with the 
planetesimal disk. Predictions of the model which 
may prove its validity are the orbital motion of 
Vega's clumpy pattern (1.1° yr^^), the location 
of the planet (8'.'3 from the star, 75° in longitude 
in front of the orbital motion of the NE clump), 
and the high resolution structure of the clumps 
including the presence of a fainter third clump on 
the opposite side of the star from the planet from 
planetesimals in the 4:3 and 5:3 resonances. 

While the mass and migration rate of Vega's 
perturbing planet are not yet fully constrained, the 
migration parameters derived for Neptune in the 
solar system (Hahn & Malhotra 1999) are close to 
the small region of parameter space which could 
have caused Vega's structure. There is also an 
intriguing suggestion that a planet's outward mi- 
gration requires the presence of massive planets 
interior to its orbit. Thus it seems possible that 
Vega's planetary system is very similar to our own, 
not only in the presence and mass of the planets 
in its system, but in that system's early evolution. 
It is also possible that application of this model to 
the other clumpy debris disks may show these to 
harbor solar-like planetary systems. The weight of 
such conclusions is damped only by the limitations 
of the model. These will be addressed in future 
studies, and include questions about the extent the 
dust distribution follows that of the parent plan- 
etesimals, and how the conclusions are affected if 
the planetesimal disk is initially dynamically hot 
(e > 0.01). 

I am very grateful to Wayne Holland for pro- 
viding the published SCUBA observations Vega. 
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A. Adiabaticity 



As a result of the fit to the hbration of 4) performed in §4.3, tlie hbration period for each planetesimal, 
t(ji was also determined for all resonant planetesimals in each run, as was the mean hbration period of the 
ensemble, {t^). The results are plotted in Figure 17. Fits to these were performed of the form {{t^) /tper )y/JI = 
A/e^ + C9/fi, where tper is the orbital period of the planetesimal and B = except for the 2:1 resonance. 
These showed that 



{h2:i)/tper = 78.9M-°■^e-o■^ (Al) 

iUsJ/tper = 317/x-°-5 - 86.40Ai^i ■^ (A2) 

iUs.J/tper = 412m-°-' - 7.70M"' ■^ (A3) 

m. J/tper = 307/x-°-^ - 5.20^"'-', (A4) 



and these fits are also shown on Figure 17. So, the hbration period for the 2:1 resonance is not affected 
by the planet's migration rate, but it does decrease during the migration as the planetesimal's eccentricity 
increases. However, for the 5:3, 3:2, and 4:3 resonances, the hbration period is constant throughout the 
migration, but does depend on the /x and 9 of the migration in that faster migrations (and larger planet 
masses) result in smaller libration periods. 

The libration width for each planetesimal, Aa = Umax — o-min, was also measured, as was the mean 
libration width (Aa). Since the runs for the 3:2 resonance showed that Aa is not just dependent on fj, and 6, 
but also on which was not varied for the other resonances, I only attempted to characterize the libration 
width for the 3:2 resonance (see also Figure 18): 

(Aa3:2) /a = 0.0023^/i-°-^^M-°-^ (A5) 

The mean libration widths for all resonances were used, however, to question the adiabaticity of this libration. 
Adiabaticity is defined as when the libration period, t^, is much shorter than the time it takes for the 
resonance to cross the libration width, i.e., tcj, <C Aa/ctpi. I found that in the runs shown here, Nm, = 
{Aa/api)/tff, was in the ranges: 1.5 — 15 for the 3:2 resonance; 2.5 — 3.9 for the 4:3 resonance; 10 — 400 for 
the 5:3 resonance; and 3.6 — 25 for the 2:1 resonance. In other words, all the migrations considered in this 
paper satisfy adiabatic invariance (Henrard 1982), though many are close to this limit. 
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Fig. 1. — Capture probabilities for the 3:2 reso- 
nance for planctcsimals initially orbiting at 60 AU 
from a 2.5Mq star plotted for different mass plan- 
ets that are migrating at different rates (a) (see 
§3.1.1). The solid lines show fits to these prob- 
abilities using the function given in equation (3). 
Parameters derived from these fits are shown in 

(b) and (c). (b) shows the migration rate re- 
quired for a planet to capture a given fraction of 
the planetesimals in its 3:2 resonance, while 7 in 

(c) defines how fast the capture probability drops 
with migration rate for a given planet mass. 
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Fig. 2. — Capture probabilities for the 3:2 reso- 
nance of a IOM0 planet for planetesimals initially 
orbiting at different distances from a 2.5Mq star 
(a) (sec §3.1.2). The solid lines show fits to these 
probabilities using the function given in equation 
(3). Parameters derived from these fits are shown 
in (b) and (c). 
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Fig. 3. — Capture probabilities for the 3:2 reso- 
nance of a IOM0 planet for planetesimals initially 
orbiting at 60 AU from stars of different mass (a) 
(see §3.1.3). The solid lines show fits to these prob- 
abilities using the function given in equation (3). 
Parameters derived from these fits are shown in 
(b) and (c). 
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Fig. 4. — Capture probabilities for planetesimals 
initially 30 AU from a 2.5Mq star for trapping 
into the (a) 2:1, (b) 4:3, and (c) 5:3 resonances. 
The solid lines show fits to these probabilities for 
each planet mass using equation (3). The dotted 
lines show the fit to the capture probabilities using 
equation (8). 
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Fig. 5. — Summary of capture probabilities for the 
4:3, 3:2, 2:1, and 5:3 resonances for migrations de- 
fined by the parameters /i and 6 (eqs. [6]-[8] with 
parameters from Table 1). The dotted, solid, and 
dashed lines indicate migrations for which trap- 
ping probabilities arc 10, 50, and 90% respec- 
tively. The dash-triple dot line indicates migra- 
tions for which the 2:l(u) resonance has twice as 
many members as the 2:1(1) resonance (see §4.3.4). 



Fig. 6. — Paths of resonant orbits in the frame co- 
rotating with the mean motion of the planet. The 
planet is marked by a cross in these figures and it 
is stationary in this reference frame, since its orbit 
is circular. Crosses show the location of a planctesi- 
mal in the different resonances (2:1, 5:3, 3:2, and 4:3) 
at intervals of l/24th of the planet's orbital period 
(i.e., when the planet has moved by 15° in the inertial 
frame). The planet moves anticlockwise in the iner- 
tial frame, and the planctesimals move anticlockwise 
in this reference frame. Sufficient crosses are marked 
to show how this pattern repeats itself, which is ev- 
ery p+q orbits of the planet. The three plots for each 
resonance show the paths for planetesimal eccentric- 
ities of 0.1, 0.2, and 0.3. Naturally, these plots show 
just one value of the planetesimal's pericenter, which 
always occurs at the same location in inertial space, 
but occurs p times in the rotating frame plots (at 
the innermost point of the loops). Resonant orbits 
with different pericenters, which can be specified by 
the resonant angle (j) (eq. [11]), would show the same 
pattern, but rotated on this figure by an angle (j)/p. 
The three plots for the 2:1 resonance have </> = 257, 
275, and 281° (for increasing eccentricity), while the 
remainder have (j) = 180° (see §4.3.3). 
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Fig. 7. Offset of the mean libration centers, 
{(j)m)i from 180° for planetesimals trapped in the 
(a) 5:3, (b) 3;2, and (c) 4:3 resonances as a re- 
sult of planet migration defined by the parameters 
H and 9 (eqs. [6] and [7]). The solid lines show 
the fits to these offsets for each resonance given in 
equations (23)- (25). 



Fig. 8. — Variation of the mean hbration center, 
{(j)m)i as a result of planet migration for planetes- 
imals trapped in the 2:l(u) and 2:1(1) resonances. 
The solid line shows the fit to this variation given 
in equation (26). The analytical solution given in 
equation (19) is shown with the dotted line. 
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Fig. 9. — Captme probabilities for the 2:1(1) res- 
onance, ^2:1(0' for planetesimals initially 30 AU 
from a 2.5M0 star. The parameterised fits to these 
probabilities given in equations (27) and (28) for 
the four planet masses shown in this plot (i.e., 
/ti = 4, 12, 40, and 120) are shown with dashed 
lines. 
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Fig. 10. — Mean libration amplitudes, (A0), of 
planetesimals captured in the (a) 2:1, (b) 5:3, (c) 
3:2, and (d) 4:3 resonances for migrations defined 
by the parameters n and 9 (eqs. [6] and [7]). The 
soHd Unes show the fits to these libration ampli- 
tudes given in equations (33)- (36). 



10* 


















A 










10^ 




















C 


















Ei 




lo" 


















^--"^ Di 










10 








Eii 








Dli 












10"** 















0.1 1.0 10.0 100.0 1000.0 

Fig. 11.- Definition of the migration zones A-Eii 
discussed in the text and Table 2. The lines are 
the same as those in plotted in Figure 5. 




0.1 1.0 10.0 100.0 1000.0 

Fig. 12. — Rotation (clockwise) of the resonant 
pattern shown in Figure 6 due to the migration of 
the planet characterized in Figure 7. This is shown 
for migrations resulting in trapping probabilities 
of 10, 50, and 90%. 
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Fig. 13. — Libration amplitudes for migrations re- 
sulting in trapping probabilities of 10, 50, and 90% 
for: (a) the 5:3, 3:2, and 4:3 resonances; (b) the 
2:1 resonance when e = 0.03 and 0.3. 
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Fig. 14. — 850 /im images of Vega, (a) Obser- 
vation taken using SCUBA at the JCMT (Hol- 
land et al. 1998). The contours start at 3.8 mJy 
beam~^ and increase at 1.9 mJy beam~^ intervals. 
The beam size has a 14" FWHM, and additional 
gaussian smoothing of 7" FWHM has been ap- 
plied. Emission from the stellar photosphere of 
'--^ 5.9 mJy has not been subtracted from the im- 
age, (b) Simulated model image of dust grains 
created in the destruction of planetesimals shown 
in Figure 15. The planet is shown at the location 
of the diamond-plus and its direction of motion 
is also shown. Appropriate color table, contours, 
smoothing and stellar photosphere have been in- 
cluded to allow a direct comparison with (a). 
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Fig. 15. — Model of the origin of Vega's structure 
as a result of the migration of a Neptune mass 
planet from 40-65 AU over 56 Myr. (a) The ini- 
tial and final orbital distribution (eccentricity vs 
semimajor axis) of planetesimals in the disk. For 
clarity only the parameters of 200 planetesimals 
(the asterisks) are shown in this plot. The planet 
is located at the diamond-plus, the dotted lines 
indicate the location of the planet's resonances, 
and the dashed lines indicate the chaotic reso- 
nance overlap region, (b) Image of the number 
density distribution of planetesimals in the disk 
at the end of the migration. The planet is located 
at the diamond-plus. 
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Fig. 16. — Alternative model explaining the 850 
/im image of Vega's disk in which the planet caus- 
ing the structure orbits the star clockwise (see Fig. 
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Fig. 17. — Mean libration periods, {t^)/tper, of 
planctcsimals captured in the (a) 2:1, (b) 5:3, (c) 
3:2, and (d) 4:3 resonances for migrations defined 
by the parameters n and 6 (eqs. [6] and [7]). The 
sohd lines show the fits to these libration periods 
given in equations (Al)-(A4). 




Fig. 18. — Mean libration width, (Aa)/a, of plan- 
etesimals captured in the 3:2 resonance for migra- 
tions defined by the parameters /U and 6 (eqs. [6] 
and [7]) for stellar masses M*. The solid line shows 
the fit to these libration widths given in equation 
(A5). 
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Table 1 

Coefficients that determine capture probabilities from equation (8) for different 

RESONANCES (§3). AlSO GIVEN ARE THE SEMIMAJOR AXES OF THE RESONANCES, Ur, WITH RESPECT TO 
THAT OF THE PERTURBING PLANET, Upi, AND THE AVERAGE ERRORS IN THE CAPTURE PROBABILITIES 

USING THESE MODELS, Pgrr- 



Resonance 


Or/ (tpl 


X 


Y 


u 


V 


Perr 


4:3 


1.21 


0.23 ± 0.02 


5.6 ± 2.0 


1.42 ± 0.01 


0.29 ± 0.1 


0.05 


3:2 


1.31 


0.37 ± 0.02 


5.4 ± 2.0 


1.37 ±0.01 


0.38 ± 0.1 


0.04 


5:3 


1.41 


210 ± 20 


1.0 ±0.2 


1.84 ±0.02 


0.20 ±0.08 


0.04 


2:1 


1.59 


5.8 ±0.2 


4.3 ±2.0 


1.40 ± 0.02 


0.27 ±0.2 


0.025 
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Table 2 

Which resonances planetesimals initially at semimajor axes of end up in given the 

INITIAL semimajor AXIS OF THE PLANET ttpi. FOR THE DIFFERENT MIGRATION SCENARIOS A-Ell SHOWN 

IN Figure 11. Also given are the longitudes of the concentrations of planetesimals in the 

DIFFERENT RESONANCES RELATIVE TO THE PLANET, A - Apj. 



Resonance 


2:l(u) 


2:1(1) 


5:3 


3:2 




4:3 


~ ^pl 
A 


-(107 - 79°) 


(107 - 79°) 


±60°, 180° 


±90° 




±60°, 180° 


B 
C 
Di 
Dii 
Ei 
Eii 


ai > 1.59api. 
50% of ai > 1.59api^ 

ai > 1.59opi^ 
50% of ai > 1.59api^ 


50% of ai > 1.59ap;. 
50% of ai > 1.59api, 


ai = (1.41 - 1.59)api. 
tti = (1.41 - 1.59)api, 


ai > 1.31ap; . 
ai = (1.31 - 1.59)ap,. 
ai = (1.31 - 1.59)op!^ 
ai = (1.31 - 1.41)api. 
ai = (1.31 - 1.41)ap,. 


ai 
ai 
ai 
ai 
ai 


ai > 1.21ap/. 
= (1.21 - 1.31)apii 
= (1.21 - 1.31)ap,. 
= (1.21 - 1.31)Opi. 
= (1.21 - 1.31)api. 
= (1.21 - 1.31)op, 
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